#include "stdafx.h"

#define FLOAT128


#include "softfloatx80.h"
#include "softfloat-round-pack.h"

static const floatx80 floatx80_negone  = packFloatx80(1, 0x3fff, IA32_CONST64(0x8000000000000000LL));
static const floatx80 floatx80_neghalf = packFloatx80(1, 0x3ffe, IA32_CONST64(0x8000000000000000LL));
static const float128 float128_ln2     = 
    packFloat128(IA32_CONST64(0x3ffe62e42fefa39eLL), IA32_CONST64(0xf35793c7673007e6LL));

#define LN2_SIG        IA32_CONST64(0xb17217f7d1cf79acLL)

#define EXP_ARR_SIZE 15

static float128 exp_arr[EXP_ARR_SIZE] =
{
    PACK_FLOAT_128(0x3fff000000000000LL, 0x0000000000000000LL), /*  1 */
    PACK_FLOAT_128(0x3ffe000000000000LL, 0x0000000000000000LL), /*  2 */
    PACK_FLOAT_128(0x3ffc555555555555LL, 0x5555555555555555LL), /*  3 */
    PACK_FLOAT_128(0x3ffa555555555555LL, 0x5555555555555555LL), /*  4 */
    PACK_FLOAT_128(0x3ff8111111111111LL, 0x1111111111111111LL), /*  5 */
    PACK_FLOAT_128(0x3ff56c16c16c16c1LL, 0x6c16c16c16c16c17LL), /*  6 */
    PACK_FLOAT_128(0x3ff2a01a01a01a01LL, 0xa01a01a01a01a01aLL), /*  7 */
    PACK_FLOAT_128(0x3fefa01a01a01a01LL, 0xa01a01a01a01a01aLL), /*  8 */
    PACK_FLOAT_128(0x3fec71de3a556c73LL, 0x38faac1c88e50017LL), /*  9 */
    PACK_FLOAT_128(0x3fe927e4fb7789f5LL, 0xc72ef016d3ea6679LL), /* 10 */
    PACK_FLOAT_128(0x3fe5ae64567f544eLL, 0x38fe747e4b837dc7LL), /* 11 */
    PACK_FLOAT_128(0x3fe21eed8eff8d89LL, 0x7b544da987acfe85LL), /* 12 */
    PACK_FLOAT_128(0x3fde6124613a86d0LL, 0x97ca38331d23af68LL), /* 13 */
    PACK_FLOAT_128(0x3fda93974a8c07c9LL, 0xd20badf145dfa3e5LL), /* 14 */
    PACK_FLOAT_128(0x3fd6ae7f3e733b81LL, 0xf11d8656b0ee8cb0LL)  /* 15 */
};

extern float128 EvalPoly(float128 x, float128 *arr, unsigned n, float_status_t &status);

/* required -1 < x < 1 */
static float128 poly_exp(float128 x, float_status_t &status)
{
/*
    //               2     3     4     5     6     7     8     9
    //  x           x     x     x     x     x     x     x     x
    // e - 1 ~ x + --- + --- + --- + --- + --- + --- + --- + --- + ...
    //              2!    3!    4!    5!    6!    7!    8!    9!
    //
    //                     2     3     4     5     6     7     8
    //              x     x     x     x     x     x     x     x
    //   = x [ 1 + --- + --- + --- + --- + --- + --- + --- + --- + ... ]
    //              2!    3!    4!    5!    6!    7!    8!    9!
    //
    //           8                          8
    //          --       2k                --        2k+1
    //   p(x) = >  C  * x           q(x) = >  C   * x
    //          --  2k                     --  2k+1
    //          k=0                        k=0
    //
    //    x                     
    //   e  - 1 ~ x * [ p(x) + x * q(x) ]
    //
*/
    float128 t = EvalPoly(x, exp_arr, EXP_ARR_SIZE, status);
    return float128_mul(t, x, status);
}

// =================================================
//                                  x
// FX2P1                   Compute 2  - 1
// =================================================

//
// Uses the following identities:
//
// 1. ----------------------------------------------------------
//      x    x*ln(2)
//     2  = e
//
// 2. ----------------------------------------------------------
//                      2     3     4     5           n
//      x        x     x     x     x     x           x
//     e  = 1 + --- + --- + --- + --- + --- + ... + --- + ...
//               1!    2!    3!    4!    5!          n!
//

floatx80 f2xm1(floatx80 a, float_status_t &status)
{
    Bit64u zSig0, zSig1;

    // handle unsupported extended double-precision floating encodings
    if (floatx80_is_unsupported(a)) 
    {
        float_raise(status, float_flag_invalid);
        return floatx80_default_nan;
    }

    Bit64u aSig = extractFloatx80Frac(a);
    Bit32s aExp = extractFloatx80Exp(a);
    int aSign = extractFloatx80Sign(a);
     
    if (aExp == 0x7FFF) {
        if ((Bit64u) (aSig<<1))
            return propagateFloatx80NaN(a, status);

        return (aSign) ? floatx80_negone : a;
    }

    if (aExp == 0) {
        if (aSig == 0) return a;
        float_raise(status, float_flag_denormal | float_flag_inexact);
        normalizeFloatx80Subnormal(aSig, &aExp, &aSig);

    tiny_argument:
        mul64To128(aSig, LN2_SIG, &zSig0, &zSig1);
        if (0 < (Bit64s) zSig0) {
            shortShift128Left(zSig0, zSig1, 1, &zSig0, &zSig1);
            --aExp;
        }
        return
            roundAndPackFloatx80(80, aSign, aExp, zSig0, zSig1, status);
    }

    float_raise(status, float_flag_inexact);

    if (aExp < 0x3FFF)
    {
        if (aExp < EXP_BIAS-68)
            goto tiny_argument;

        /* ******************************** */
        /* using float128 for approximation */
        /* ******************************** */

        float128 x = floatx80_to_float128(a, status);
        x = float128_mul(x, float128_ln2, status);
        x = poly_exp(x, status);
        return float128_to_floatx80(x, status);
    }
    else 
    {
        if ((a.exp == 0xBFFF) && (! (aSig<<1)))
           return floatx80_neghalf;

        return a;
    }
}
